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Abstract. We show that the continuous wavelet transform can provide a unique de- 
composition of a timeseries in to 'signal-like' and 'noise-like' components: From the overall 
wavelet spectrum two mutually independent skeleton spectra can be extracted, allowing 
the separate detection and monitoring in even non-stationary timeseries of the evolution 
of (a) both stable but also transient, evolving periodicities, such as the output of low 
dimensional dynamical systems and (b) scale-invariant structures, such as discontinuities, 
self-similar structures or noise. The idea of the method is to keep from the overall wavelet 
expansion of the timeseries only the wavelet components of largest amplitude at any given 
time or scale, thus obtaining the instantly maximal and scale maximal wavelet skeleton 
spectrum, respectively. The scale maximal spectrum was previously proposed for studying 
possible multifractal scaling properties of the timeseries (e.g. Arneodo et al, 1988). The 
here proposed instantly maximal spectrum exhibits clearer spectral peaks and reduced 
noise, as compared to the overall wavelet spectrum. An indicative application to the 
monthly-averaged sunspot index reveals, apart from the well-known 11-year periodicity, 3 
of its harmonics, the 2-year periodicity (quasi-biennial oscillation, QBO) and several more 
(some of which detected previously in various solar, earth-solar connection and climate in- 
dices), here proposed being just harmonics of the QBO, in all supporting the double-cycle 
solar magnetic dynamo model (Benevolenskaya, 1998, 2000). The scale maximal spec- 
trum reveals the presence of 1// fluctuations with timescales up to 1 year in the sunspot 
number, indicating that the solar magnetic configurations involved in the transient so- 
lar activity phenomena with those characteristic timescales are in a self-organized-critical 
state (SOC), as previously proposed for the solar flare occurence (Lu & Hamilton, 1991). 

1 Introduction 

The widely used Fourier transform, although useful for stationary signals of sim- 
ple dynamics consisted of a linear superposition of few independent, strong, non- 
evolving periodicities, has severe drawbacks for analyzing signals of the following two 
important categories, often encountered as output of physically interesting complex 
systems: 

a. Signals that include transient or variable periodicities. The Fourier transform, 
since integrating over the whole time domain, does not allow monitoring of amplitude 
or frequency evolution in periodicities. A limited solution is to divide the signal in 
segments or else 'windows' (a method called short-time, Gabor or windowed Fourier 
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transform), however introducing the arbitrary, fixed window width that imposes a 
low frequency cut-off near which edge effects occur, thus artificially deforming the 
estimated spectrum. 

b. Signals that significantly depart from stationarity, consisted of intermittent 
'activity bursts' (e.g. energy releases), or hierarchies of localized structures in space 
or time (e.g. in turbulent fluids reside eddies, discontinuities, filaments, sheets, 
shocks etc.). For each localized structure a wide frequency range of strongly phase 
coherent (i.e. 'synchronized' in time) Fourier sine components are required to repro- 
duce it, adding near and canceling away from the structure. Hence, each structure 
spreads and become undetectable over the whole Fourier spectrum. Moreover, such 
structures often exhibit 'self-similarity', or else 'scale invariance', i.e. consists of 
hierarchies of mutually similar structures, emerging through a cascade mechanism, 
and 'roughness' over a wide scale range (well modeled by fractals and multi-fractals) 
that is not built in the sine waves of Fourier transform. 

The projection of signals with continuous or discrete wavelet analyzing func- 
tions, that are by definition self-similar and well localized in both time (or space) 
and frequency (or wavenumber) domains, provides efficient descriptions, in terms 
of sparseness, of non-stationary, transient, or 'bursty' records. It is proven a useful 
tool, used in an increasingly large variety of applications. Here we mention some of 
them indicatively according to the type of use: detection of periodic signals in noisy 
timeseries (e.g. Otazu et al, 2002), monitoring of period variations (e.g. Frick et 
al, 1997, Fligge et al, 1999), unevenly sampled period analysis (e.g. Foster, 1996), 
derivation of multifractal properties (e.g. Arneodo et al, 1988, Argoul et al, 1989, 
Muzy et al, 1991, 1993, Arneodo et al, 1995, 1998), random multifractal synthe- 
sis (e.g. Benzi etal, 1993), localized structure identification (e.g. Lucek & Balogh, 
1997, Roux et al, 1999, Mouri et al, 1999), polarization analysis (e.g. Baumjohann 
et al, 1999), denoising (e.g. Fligge & Solanki, 1997, Komm et al, 1999), signal- 
to-noise ratio enhancement (e.g. Zhang & Paulson, 1997), information compression 
(e.g. Muhlmann & Hanslmeier, 1996), multi-resolution image decomposition (e.g. 
Mallat, 1989, Pantin & Starck, 1996), studies of random walks (Arneodo et al, 
1996), flow structure analysis (Haynes & Norton, 1993), population classification 
(Bendjoya et al, 1991), clustering detection (e.g. Bijaoui et al, 1993, Girardi et al, 
1997, Lima Neto et al, 1997), fast data query (e.g. Chakrabarti et al, 2001), fluid 
turbulence simulations (e.g. Schneider et al, 1997) etc. 

The wavelet transform is now used extensively for the analysis of oscillations 
and periodicities in various solar structural and activity features, as the sunspot 
number (Fligge et al, 1999, Mordvidov & Kuklin, 1999, Sello, 2000, Rigozo et 
al, 2001, Nayar et al, 2002), sunspot groups (Ballester & Oliver, 1999), active 
regions (Ireland et al, 1999, O'Shea et al, 2001, 2002), coronal holes (Marsh et al, 
2002), bursts (Schwarz et al, 1998, Meszarosova et al, 1999, Gimenez et al, 2001), 
flares (Aschwanden et al, 1998), polar plumes (Banerjee et al, 2000), coronal loops 
(De Moortel et al, 2000, 2002a, b), photospheric flows (Lawrence et al, 2001), 
chromosphere (Bocchialini & Baudin, 1995, Frick et al, 1997), transition region 
(Fludra, 2001), differential rotation (Hempelmann & Donahue, 1997, Soon et al, 
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1999, Hempelmann, 2002) etc. 

In this article we demonstrate that two exactly defined, mutually independent, 
complementary wavelet skeleton spectra, extracted from the total continuous wavelet 
spectrum, allow to discriminate the components of (stable or transient) periodici- 
ties and of hierarchies of discontinuities in timeseries. An indicative application is 
made to sunspot number, since (a) is known to include both periodic and noise-like 
components, and (b) its study may contribute in understanding the features and 
mechanism of the solar activity. Section 2 includes a brief review of the continuous 
wavelet transform for the unfamiliar reader (for extended reviews see e.g. Hunt et 
al, 1993, Kaiser, 1994, Louis et al., 1997, Torrence & Compo 1998). The two wavelet 
skeleton spectra are introduced in Section 3. The application to the sunspot index 
is presented in Section 4. Conclusions and discussion are summarized in Section 5. 

2 A brief review of the wavelet transform 
2.1 Linear functional projection 

One way to quantify the degree of 'similarity' between two, generally complex func- 
tion x(t), y(t) of a variable t G [0, T] (e.g. time or distance) is by the amplitude (or 
else cross-correlation) of their functional projection over the domain of t: 



with y* being the complex conjugate of y. The factor 1/T makes the amplitude a 
independent of T and finite in the limit T — > +oo. The two functions are called 
'orthogonal' if a = 0. 

In analogy to the projection of a vector to a coordinate system in physical space, 
a given function x(t) (e.g. the timeseries of a measured observable) thus can be 
analyzed to a family of functions y(t,Pi), of n characteristic parameters pi, i = l...n: 



The functional analysis (2) is linear in the sense that any x(t) can be expanded to 
a linear superposition of a set of functions y(t,Pi) that are mutually orthogonal and 
integrable. 

Note that the continuous functional transform (2) is generally redundant since it 
spreads the information content of x(t) from the one-dimensional time axis, to the 
n-dimensional space of the parameters Pi. The minimization of information spread 
in the parameter space pi by appropriate selection of the functional basis and the 
development of efficient redundancy reduction schemes is of central importance for 
spectral analysis, functional approximation and compression applications. 




(1) 




(2) 
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2.2 Appropriate selection of analyzing functions 

Although the linear functional projection (2) is mathematically valid with any set of 
functions, an appropriate projection should aim in selecting such functions to match 
as much as possible (i.e. being nearly orthogonal to most of the functions so that 
the representation is sparse) the analyzed function, which often is or consists of the 
following: 

(a) Analog signal, as the independent variables of a system or model of relatively 
simple dynamics. The usually smoothly varying variables of systems of low dimen- 
sionality (i.e. small number of independent variables) exhibit oscillatory behavior, 
thus can often be sparsely expanded to a linear superposition of a few periodic 
functions, or, more generally, transient, quasi-periodic components. Hence, it is 
desirable that the analyzing functions have well-defined, single-peaked, narrow fre- 
quency spectral content. 

(b) Noise, as this due to observational 'random' fluctuations (often additive) 
of high-frequency content, e.g. of flat spectrum in uncorrelated (white) or more 
generally of wide range power-law decaying spectrum in auto-correlated ('colored') 
noise. Noise-like contributions may also caused by fast-varying components of a 
deterministic complex system of large number of independent variables (such as 
macroscopic systems) or the interaction of a low-dimensional open system with an 
environment of high dimensionality. The observables of high-dimensional systems, 
often involving strongly turbulent fluids, are characterized by intermittency or 'sin- 
gularities' due to coherent structures or intermittent energy bursts and, exact or 
weaker, self- similarity between different scales (often termed scale-invariance) , well 
modeled by fractals and multi-fractals. Although such systems are characterized 
by strongly non-linear physical mechanisms, for an appropriate linear projection 
of such observables, the analyzing functions should at least posses the element of 
self-similarity. 

2.3 Fourier transform 

The widely used Fourier transform can be interpreted as mapping, or complete 
'rotation' (in analogy to the vector rotation in physical space) from the time domain 
to the frequency domain, by means of the complex periodic plain wave functions: 



thus represents an estimate of the amplitude of periodicities of frequency /, that 
occurred during the time interval T. The, often studied, Fourier power spectral 
density is defined as follows: 




(3) 




(4) 



SIGNAL-NOISE DECOMPOSITION: APPLICATION TO SUNSPOTS 



5 



P(f)=a 2 (f) (5) 

A useful property of the Fourier transforms is the power theorem: any linear func- 
tional projection (2) can be calculated via Fourier transforms in the frequency do- 
main (e.g. Hunt et al, 1993, p. 6): 

rT r+oo 

a(pi) = / x(t) y*(t, Pi ) dt = / x(f) y*(f, Pi ) df (6) 

JO J-oo 

where x, y the Fourier transforms of x, y respectively (see (4)). Hence, the functional 
projection (2) can be treated as a linear filter acting on x. 

2.4 Time- frequency analysis 

Note that each Fourier basis function y is totally localized in frequency but nowhere 
in time, since having infinite duration (i.e. y in (3) do not vanish for t — > ±oo). As a 
consequence, the Fourier transform does not allow monitoring of transient (periodic 
or aperiodic) components in the signal, since they are spread over the whole Fourier 
spectrum. 

In the framework of the transform (2) a solution to this problem is to introduce 
analysing functions of well-defined frequency but are also localized in time i.e. of 
significant amplitude for a finite time duration. 

Perfect localization in both time and frequency in however inhibited by the 
uncertainty principle (e.g. Kaiser, 1994, p. 50): (2) keeps the information of x(t) 
confined in cells, within which there is total uncertainty about the information of 
x(t), and the size of these cells cannot be infinitesimally small, having a lower limit 
depending on the properties (shape, rate of decay etc.) of the selected analysing 
functions. 

Any time-frequency transform represents a mapping of x(t) on the t — f plane. 
Given a set of (generally complex- valued) functions y(t,t',f), of characteristic fre- 
quency /, localized in the vicinity of time t', the linear functional projection (2) 
is: 

a(t,f) = ± J\{t')y*{t,t',f)dt' (7) 

The graph of the power a 2 (t, f) on the t — f plane is oftenly called scalogram. The 
total spectrum: 

p(f)= [ T a 2 (t,f)dt (8) 
J o 

is expected to generally exhibit similar features and shape (e.g. spectral peaks, 
slope) as the corresponding Fourier spectrum (equ. (5)). 
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2.5 Windowed Fourier transform 

In the windowed Fourier transform (Gabor, 1946) the Fourier plain wave functions 
are modulated by a 'support' (often called 'window' function) of significant value 
only in the vicinity a of time t': 

y^t>J) = e^ t <p( t —) (9) 

a 

Its Fourier transform has spectral width of the order of 1/a. Therefore, the t — f 
plane is resolved in cells of fixed shape and area at least a by 1/a. 
A commonly used support function is a Gaussian of width a: 

,,t-t\ (t-t'f , . 

0( ) = e~^^~ (10) 

a 

The resulting analysing function (9) has Fourier transform: 

y(fj , )=*e= i£ ^ 1 (11) 

The Gaussian support can be proven to give optimal localization in the sense that 
among all functions it covers smallest area permitted by the uncertainty principle 
in the time-frequency plane (Louis et al., 1997). Other window functions are also 
used, aimed at reducing the leakage of power of spectral peaks to nearby frequencies 
(e.g. Press et al, 1992). 

However the windowed Fourier transform always introduces the arbitrary dura- 
tion of the support function (a in equ. (10)) which itself does not depend on / and 
must be selected appropriately before the analysis. Also, the analyzing functions y 
are not self-similar since include increasingly many oscillations within their support 
with increasing frequency. 

2.6 Continuous wavelet transform 

The continuous wavelet transform resolves the information on the t — f plane in cells 
of variable shape depending on frequency, with size of the order of 1/f by /. With 
this partition of the t — f plane the arbitrary selection of the windowed Fourier trans- 
form is avoided and the analyzing functions (wavelets) are exactly self-similar (or else 
scale-invariant) since each includes a constant number of oscillations independent of 
frequency. The wavelet transform at lower frequencies provides better resolution in 
frequency but worse in time (since the analysing functions are wider) and at high 
frequencies better time localization (since the analysing functions are narrower) but 
more uncertainty in frequency, in all acting as a 'mathematical microscope'. 

While many wavelet families are so far proposed (also depending on the nature 
of application), of the most commonly used in physical applications is the Mor- 
let wavelet (Morlet et al, 1982), being the generalization of the windowed Fourier 
transform (see equ. (9)): 
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(12) 



-/ 2 (t-t') : 

#/(f-0) = e— ^ 



(13) 



Its Fourier transform is also Gaussian of width /: 



1 (/-/') 

y(fJ') = 7 e "^ 



(14) 



Among all possible wavelet analysing functions, the Morlet wavelet due to its Gaus- 
sian support inherits optimality as regarding the uncertainty principle (Louis et al, 



3 The two wavelet skeleton spectra 
3.1 Signal-noise decomposition: 

Successful decomposition of the timeseries into 'signal' and noise-like component is 
of fundamental importance in various applications, for the separate study of the 
signal and noise properties (e.g. locating spectral peaks and deriving scaling laws), 
denoising, compression etc. The wavelet transform is promising for that decompo- 
sition, since allowing monitoring of the evolution of time-localized wavepackets of 
specific frequency, hence also of stable or transient periodicities and noise-like spikes. 

The overall wavelet scalogram generally contains a mixture of the signal and 
noise characteristics. Imposing some thresholding logic (e.g. hard or soft level, 
trend etc) to the wavelet components (i.e. a type of wavelet band-pass filter) is 
expected to have limited success in the case of strong, auto-correlated noise and/or 
weak, transient periodicities. The deeper reason that such filters have limited success 
in those cases is that the selection of the wavelet components that are attributed 
to the signal or the noise is to some extend arbitrary, not adaptively taking into 
account the information content in the original timeseries. 

As commented in section (2.1) for any continuous functional projection, the con- 
tinuous wavelet transform is always redundant in both time and frequency, since 
mapping the one dimensional information in x(t) on the (two dimensional) t — f 
plane. It contains correlations between the wavelet components that do not exist 
in the timeseries but is significant at those (small) regions and scales where the 
wavelet functions cross-correlation is also significant. As a consequence, the higher 
the achieved resolution (i.e. sampling in frequency), the higher the smoothness of 
the resulting spectrogram due to the increased cross correlations between wavelet 
amplitudes at nearby frequencies. The presence of noise-like components, e.g. due 
to observational noise, contaminates ('blurs') the scalogram, by making more un- 
correlated wavelet amplitudes of nearby times and introducing ridges across wide 
ranges of frequencies. 



1997). 
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In contrast to the continuous, the discrete wavelet transform (the review of which 
is outside the scope of this article) leads to uncorrelated wavelet amplitudes if the 
discrete wavelet analysing functions are orthogonal. However the lack or redundancy 
is not an advantage for the signal-noise decomposition, since there is no profound cri- 
terion for attributing the wavelet components to signal or to noise. Some threshold- 
ing criterion can be imposed, as commonly used in signal compression applications, 
again having a degree of arbitrariness. 

From the above discussion emerges a plan for an efficient signal-noise decompo- 
sition using the continuous wavelet transform, consisted of the following steps: 

(a) Continuous wavelet transform of the timeseries with sufficiently high fre- 
quency resolution (sampling), thus producing a highly redundant wavelet scalogram. 

(b) A non-parametric, information dependent selection criterion that disposes 
effectively the redundancy, in a way that selects separately those wavelet components 
that can be with high probability attributed to signal and to those to noise. 



3.2 Scale maximal wavelet skeleton spectrum 

A transient structure or burst of finite, characteristic duration r, occurring (i.e. 
having center of mass) in the timeseries at time t is expected to cause a local increase 
of the amplitudes \a(t, f) \ of the timeseries wavelet transform in the vicinity of time 
t, at wavelet period of the order of r (or else frequency 1/r), while can extent 
across different timescales, e.g. a discontinuity (r = 0) is wavelet-transformed to a 
superposition of wavelets with the same center of mass at time t, at all timescales 
(and frequencies). 

For the study of such structures the scale maximal wavelet skeleton spectrum is 
extracted from the overall wavelet spectrum, keeping only those wavelet components 
of which the complex amplitude is locally maximum across time at any given time- 
scale, i.e. those for which: 

*M^0, d ^ m <0 (15, 
at ot z 

It is already recognized that the so defined spectrum (otherwise called wavelet trans- 
form modulus maxima) contains all the important information about the existence 
of hierarchies of discontinuities, usually appearing as continuous lines across scales 
pointing at the time of the discontinuity, hence can be used for detection of singu- 
larities (e.g. Mallat, 1999, p. 176), or even approximate reconstruction of timeseries 
(e.g. Mallat, 1999, p. 185 and references therein). It is also used for the study scal- 
ing (e.g. possible multifractal) properties in timeseries: the lines of each singularity 
is often, at least at high frequencies, consisted of wavelet amplitudes that follow 
power-law dependence on frequency across many orders of frequency range (i.e. for 
these ranges there is no characteristic scale, or else there is some scale-invariance 
in the timeseries), well modeled by multifractals: If the timeseries is a realization 
of a multifractal then the wavelet amplitudes of the locally scale-invariant skeleton 
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spectrum, at least at high frequencies, / — ■> +00, follow power-law of the form (for 
a review see e.g. Arneodo et ai, 1995): 

|a(f,/)|<x/- fc « (16) 

where h(t) is the Holder (or else local Hurst) exponent, being a local measure of 
'burstiness' in the timeseries at time t (for details see e.g. Muzy et ai, 1993). 

Power-law spectral power dependance over frequency, at least over several orders 
of frequency magnitude and especially tails at high frequencies (or small scales) are 
fairly common in all physical observables. In the special case of uni-fractals h(t) 
has the same value for all t. Of specific interest is the case of 1/ f noise (i.e. with 
h — 0), constituting the hallmark of self- organized criticality, (SOC, Bak et ai, 
1987) i.e. the meta-stable state close to a critical point of a universal class of systems 
characterized by intermittency, well modeling aspects of solar flares, earthquakes, 
geological formation, biological evolution, economy etc. (for a review see Bak, 1996). 



3.3 Instantly maximal wavelet skeleton spectrum 

A transient or stable periodic component of instant period T at time t in the time- 
series is expected to cause increase of the wavelet amplitudes \a(t, f) \ in the vicinity 
of wavelet period of the order of T (or else frequency 1/T), while may extend across 
time as far as its period T remains constant in time. 

For the detection of periodicities in timeseries here we introduce the instantly 
maximal wavelet skeleton spectrum, extracted from the overall wavelet spectrum, if 
keeping only those wavelet components that are locally of maximum amplitude at 
a given time, i.e. only those for which: 

*MM = , ^Ml <0 (17) 

df dp 

From the above discussion becomes clear that the so defined spectrum can be used 
for detecting transient periodicities of varying period, and of course those of stable 
period in the timeseries. Note that the instantly maximal skeleton spectrum excludes 
the wavelet components that instantly extend across time scales, thus it is less 
contaminated by the noise-like scale-invariant component of the timeseries. On the 
other hand, the scale maximal spectrum includes only those wavelet components that 
extend across time scales, thus the scaling properties derived from it are not distorted 
by possible transient periodicities in the timeseries due to the presence of components 
of analog signal type. In that sense the two skeleton spectra are complementary, 
permitting the separate study of the noise-like and signal-like timeseries component 
properties separately. Noise-like components of an hierarchy of discontinuities also 
affects strongly the complex phases, arg(a(t, /)), of the wavelet spectrum since the 
wavelet expansion requires 'synchronization' of wavelets (i.e. to be placed nearby in 
time). For this reason only the wavelet amplitudes are considered in (17). However, 
auto-correlated random noise has significant slowly varying components, i.e. non- 
periodic trends at various scales (and in that sense the two skeleton spectra do 
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not completely separate signal from noise), nevertheless introducing a power-law 
background without concentration of power in spectral peaks. In all, the instantly 
maximal skeleton spectrum is expected to exhibit more clearly the periodicities of 
a signal-like component of timeseries, practically excluding uncorrelated hierarchies 
of discontinuities due to noise. 

4 Application to the sunspot index 

4.1 A review on solar activity and sunspots 

The solar activity (e.g. sunspots, flares, coronal mass ejections etc) is rich, involving 
complex, not well understood magnetofluid processes, extending in a wide range of 
spatial scales and durations. The sunspot number is widely regarded as an easily 
collectable quantitive solar activity measure. Sporadic, naked-eye observations exist 
in Chinese dynastic histories since 28 BC, while systematic sunspot observations are 
kept since the discovery of telescope in 1610, while not uniformly reliable (Waldmeier 
1961, Eddy, 1977, Sonett 1983). 

The sunspot number is dominated by the 11-year periodicity (well-known Schwabe 
in 1843) connected to the 22-year cycle of the overall solar magnetic field. The 11- 
year cycle is stable, exhibiting only small variations of each cycle duration, as also 
shown using wavelet analysis of sunspot number (Fligge et al, 1999). However, the 
maximum number of sunspots in each cycle varies significantly and quite unpre- 
dictably, even nearly, but not totally, vanishing for many cycles. In all, the sunspot 
record includes enigmatic epochs of suppression, as the Maunder Minimum from 
1645 to 1715 (e.g. Ribes & Nesme-Ribes, 1993) and Sporer Minimum in the 15 </l 
century, as well as epochs of enhanced activity, as the nowadays Modern Maximum 
and Medieval Maximum in the 12 th century. The reason for the large variations of 
the cycle's amplitudes and the epochs of suppressed activity of Maunder minimum 
type is unclear, e.g. proposed being due to chaotic behaviour of the nonlinear dy- 
namo equations (e.g. Ruzmaikin, 1983, Kiiker et al, 1999) or due to stochastic 
instabilities forcing the solar dynamo leading to on-off intermittency (e.g. Hoyng, 
1988, 1993, Ossendrijver & Hoyng, 1996, Schmitt et al, 1996). 

The average shape of each sunspot cycle is systematically asymmetric, taking 
less time rising to maximum than reaching the next minimum, implying that the 
solar cycle is intrinsically non- linear (e.g. Veselovski & Tarsina, 2002). Probably the 
most significant relation between the cycle shape characteristics is the Waldmeier 
effect: cycles with larger amplitudes are more asymmetric, taking less time to reach 
maximum (e.g. Hathaway, et al, 1994, Li, 1999, Veselovski & Tarsina, 2002). That 
the 11-year oscillation is not harmonic has an important implication: it introduces 
infinite number of spectral peaks (harmonics) with periods 11/n years, n = 2,3, ... 
(e.g. Polygiannakis et al, 1996, Mursula et al., 1997). All those general features 
were successfully described by a simplified mono-parametric Van der Pol non-linear 
RLC oscillator model, further shown to accurately reproduce the observed sunspot 
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record and extended epochs of suppressed activity, given the maximum of each cycle 
(Polygiannakis et al, 1996). 

Systematic search for periodicities other than the basic 11-year in sunspot num- 
ber, solar activity indices and, consequently, in the interplanetary, geomagnetic, 
earth rotation fluctuations and climate parameters lead to numerous claims, also 
depending on the nature of the analyzed parameter (as regarding its sensitivity to 
the solar activity periodicities), the periodicity detection method and the finite size 
of the records: 

(a) Periodicities longer than 11 years, such as 180 years, 90 years ('Gleiseberg 
cycle'), 45 and 22 years (e.g. Cohen k Lintz, 1974, Sonnet 1982), perhaps subhar- 
monics (i.e. integer multiples) of the basic 11-year periodicity, and several others 
(Norderman k Trivedi, 1992), intended to describe the long-term modulation of 
the sunspot cycle amplitudes as periodic, in contrast to the chaotic or stochastic 
theoretical scenarios discussed above. 

(b) 11-year periodicity harmonics, namely of 5.5 (FYO, five-year oscillation) and 
3.7 years in sunspot number, other solar activity and solar-terrestrial connection 
parameters as the geomagnetic field and earth rotation disturbances (e.g. Currie, 
1976, Courtillot k Le Mouel, 1976, Suguira, 1980, Carta et al, 1982, Kondor, 1993, 
Djurovic k Paquet, 1996, Mursula et al, 1997, Nayar et al, 2002). 

(c) Numerous shorter periodicities of transient character, detected over one or 
few 11-year cycles, with most frequently appearing in bibliography the 2-year (QBO, 
quasi-biennial) oscillation in sunspot number (Shapiro k Ward, 1962, Apostolov, 
1985), solar radio flux (Hughes k Kesteven, 1981), solar magnetic helicity (Bao 
k Zhang, 1998), other solar activity, solar-terrestrial connection and atmospheric 
parameters and earth rotation rate (e.g. Chao, 1989, Djurovic k Paquet, 1993, 
Kane, 1997). Oftenly discussed are also periodicities of 154, 129, 103, 77, 51 days 
in solar flare occurrence and sunspot area (Riegel et al, 1984, Bogard k Bai, 1985, 
Dennis, 1985, Kile k Cliver, 1991, Bai k Surrock, 1991, Carbonel k Ballester, 1992). 

The 11-year sunspot (equivalently, 22-year solar magnetic) cycle phenomenon, 
is now believed to be due to a magneto-hydrodynamic dynamo action, periodically 
regenerating the solar magnetic field at the basis of the convection zone. The cause 
and importance of periodicities other than the 11-year remain so far unclear. The 
periodicities of 154, 129, 103, 77, 51 days were proposed to be subharmonics of a 
25.5 days fundamental periodicity of some 'clock mechanism' in the sun, modeled 
by an oblique rotator, (Bai k Surrock, 1991, Sturrock k Bai, 1992). Alternatively it 
was proposed that the 154-day periodicity and others can be related to the normal 
modes of the solar interior oscillations (Wolff, 1983, 1992). The source of the QBO 
may be connected with periodic poleward streams of magnetic flux, first discussed 
by Howard and LaBonde (1981). Based on magnetic helicity observations (Bao k 
Zhang, 1998), it was proposed a double-cycle solar magnetic dynamo model, with one 
dynamo operating at the basis of the convection zone (11-year cycle) due to radial 
shear and the other at the top of the convection zone (with 2-year periodicity) due 
to latitudal shear (Benevolenskaya, 1998, 2000). Also supporting this model, further 
analysis of the solar magnetic field data showed that the active region magnetic field 
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exhibits an intense, short-scale component varying with 2-year periodicity (Erofeev, 
2001). It was also proposed that the sunspot number can be well fitted by a super- 
position of the usual 11-year cycles and wave trains with periodicity continuously 
varying from 38 months at solar maximum, to 21 months towards solar minimum 
(Krasotkin, 2001). 

The sunspot number and all other solar activity parameters also exhibit non- 
periodic intense fluctuations that should be attributed to the episodic character 
of short-timescale solar phenomena. The solar magnetic field on the photosphere 
exhibits complex organization and sudden, intermittent energy emissions as the solar 
flares, in all well modeled as if the solar corona is in a self-organized critical state 
(SOC). The corresponding cellular automata models are in good agreement with the 
observed power-law distributions of flare released energy and occurrence rate (e.g. 
Dennis, 1985, Lu & Hamilton, 1991, Lu et al, 1993, Georgoulis & Vlahos, 1998, 
Isliker et al, 2000, 2001, Anastasiadis, 2002). The emergence of a flare drastically 
reorganizes the magnetic fields in solar active regions, hence is expected to affect 
the underlying sunspot and active region lifetimes, thus leaving an imprint of a, 
characteristic of SOC, 1/ / tail at high frequencies in all the solar-activity related 
parameters, as the here considered sunspot index. 

From the timeseries point of view the sunspot number record provides a good 
example of mixture of, both stable and transient, periodicities of significantly varying 
amplitudes and of high-frequency, noise-like components (e.g. Watari, 1996). From 
the physical point of view it promises to provide useful information about the, not 
completely understood, nature of the solar activity, as connected to the overall solar 
magnetic field regeneration mechanism and, at smaller timescales, the generation of 
active, transient solar phenomena, such as flares and coronal mass ejections. 

4.2 Wavelet spectra of sunspot index 

For this application we the record of 2048 recent monthly sunspot number timeseries 
extending from June 1831 to February 2002, compiled from the SIDC RWC Belgium 
World Data Center for the Sunspot Index catalogues. The number of points was 
selected being a power of 2 in order to avoid edge effects in the fast Fourier transforms 
used to derive the wavelet scalograms, thus also omitting the record before 1830 
which are of questionable reliability (e.g. Sonett 1983). The use of monthly-averaged 
numbers excludes effects related to the solar differential rotation, with period of 27 
days at the solar equator, and the time-life of each sunspot, persisting for one or for 
a few solar rotations. 

The wavelet transform was calculated using the Morlet wavelet (equ. (13)), 
filtering the timeseries for each wavelet scale in the frequency domain for N = 100 
total wavelet periods, sampled logarithmically between the minimum sampling time 
(1 month) and the total timeseries duration (2048 months) and then inverse Fourier 
transforming it in the time domain (power theorem, see (6)). Note that the total 
number of period samples, N, used to derive the wavelet scalogram is arbitrary, 
depending on the information content in the record (defining the complexity of the 



SIGNAL-NOISE DECOMPOSITION: APPLICATION TO SUNSPOTS 



13 



scalogram). Too large N may lead to mode splitting effects and too small to poorly 
sampled scalogram. Generally, selecting N being of the order of 1/10 of the total 
number of points in the timeseries (here 2048) gives satisfactory results. For a review 
on continuous wavelet transform and the details of the related numerical techniques 
used to derive the wavelet scalogram, see Torrence & Compo (1998). Once the 
wavelet scalogram is calculated it is parsed twice, once across time at each scale and 
once across scales for each time instance, testing successive triads of amplitudes and 
keeping only the locally maximal, thus obtaining the two skeleton scalograms. 

In figure 1(a) we present the monthly- averaged sunspot index. The well-known 
11-year periodicity and the asymmetry of each cycle are prominent. Also, the fluc- 
tuations of the index are larger near each cycle's maximum. 

In figure 1(b) the corresponding wavelet scalogram is shown. Notice the stable 
11-year periodicity, having nearly constant period with time. The spectrum is rich 
in higher-frequency components too, and the intensity of higher frequencies is larger 
at cycle's maxima, since the fluctuations are more intense there. Note that there is 
some arbitrarily in the image of the scalogram (as this applies to the representation of 
any wavelet scalogram and more generally time- frequency distribution), in the sense 
that appropriate color or grey scales must be attributed to the range of wavelet 
amplitudes, resulting in visually richer or poorer images. 

The two skeleton scalograms, extracted from the overall wavelet one, are shown 
in figures 1(c) and 1(d). In the locally scale- invariant spectrum it can be noticed that 
the lines of maxima proceed towards larger periods at the cycle's maxima, since they 
represent the center of each cycle as an overall structure. The instantly maximal 
skeleton spectrum clearly shows the 11-year periodicity as a nearly horizontal line 
and others, also quite stable. Periodicities of about 22-years are also prominent, 
probably being sub-harmonics of the 11-year one. Near the period of 1 year the 
periodicities appear having significant evolution with time, lasting for many decades 
(several cycles). The periodicities of smaller periods seem becoming increasingly 
unstable and intermittent (i.e. the lines have more and more gaps), probably due 
to their small amplitude, and the presence of self-similar noise-like fluctuations. 

The total wavelet power of the overall and the two skeleton spectra is shown in 
figure 2. The main spectral peaks are indicated with arrows, noting the correspond- 
ing period in years. On the first column of table 1 we give the period of each of 
these peaks. The corresponding period error estimates represent the width between 
successive wavelet periods, as sampled to derive the wavelet transform (equ. (13)). 

The overall wavelet spectrum is poor, exhibiting only few spectral peaks, namely 
at the 11, 5.5, 2 and 1 year, and an 1/ / tail at periods smaller than 1-year. However, 
the instantly maximal skeleton spectrum exhibits several spectral peaks, and this 
verifies its value: apart from the dominant 11-year periodicity, 3 of its harmonics, 
expected due to the 11-year cycle rise-fall asymmetry, with periods 11/n, n = 2, 3, 4 
(i.e. 5.5, 3.67 and 2.75 years) and rapidly decreasing amplitudes at shorter periods 
are also prominent. The weaker peaks of 16, 13 and 8.8 years seem related with the 
main 11-year periodicity, probably being a mode-splitting effect. The long, 28,24 
and 19 year periodicities are probably due to mode-splitting of a 22-year periodicity, 
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being the main 11-year sub-harmonic, or can be interpreted as modulation effects 
of the cycle's amplitudes (see the discussion for periods longer than 11-years in the 
previous section). 

At shorter periods, a set of 10 spectral peaks with less rapidly decreasing am- 
plitudes (nearly as 1/f with frequency) exists, with periods given in Table 1. Note 
that the QBO and the periodicities of 154, 129, 103, 77 days, i.e. 0.42, 0.35, 0.28, 
0.21 years (previously reported in various solar activity parameters, as discussed in 
the previous section), are recovered. 

As discussed in the previous section, according to the 'clock mechanism' scenario 
(Bai & Surrock, 1991, Sturrock & Bai, 1992) the 154, 129, 103, 77 day periodicities 
constitute the 6 th , 5 th , 4 th and 3 th subharmonics of a 25.5 days fundamental periodic- 
ity. However the remaining 6 peaks are not easily explained in a similar way as even 
higher subharmonics. The QBO should be the 29 th subharmonic of the 25.5 day 
periodicity, however the subharmonics in between are not observed in the spectrum. 

Alternatively, here we propose that all these periodicities are actually harmonics 
of the QBO periodicity, which then itself should be non-harmonic. In table 1 we 
present the periodicities found and also those predicted by the hypothesis of only 
two non-harmonic fundamental oscillations, of 11 and 2 years. Note that 12 out 
of 15 detected spectral peaks, with periods of 11-year and shorter can be explained 
solely with this hypothesis. However: 

(a) Two of the predicted QBO harmonics, of 0.29 and 0.22 years are not promi- 
nent in the spectrum, possibly being merged, due to noise, with the expected 0.33 
and 0.2 year peaks respectively. 

(b) The expected single spectral peak at 1-year period is wide, extending from 
about 0.9 to 1.6 years, within which 3 merged peaks of 0.97, 1.07, 1.4 period are 
detectable. Possibly this peak widening is due to evolution of ephemeral, evolving 
periodicities (as seen in fig. 1(d)), indicative of transient underlying solar activity, or, 
alternatively, the 0.97 and 1.4 year peaks are additional periodicities (independent 
of the QBO) of unknown origin. 

The scale maximal skeleton spectrum is affected only by the main 11-year pe- 
riodicity, appearing as a dominant spectral peak, while at timescales of 1-3 years, 
where the main QBO peak resides, the scale maximal spectral power is reduced. At 
periods smaller than 1 year the scale maximal spectrum exhibits a plateau of nearly 
zero-slope. This power law-behavior at high frequencies is expected from (16) for 
Hurst exponent h = 0. The noise-like component in the sunspot index thus exhibits 
power-law (scale-invariant) properties characterizing the 1/f noise, supporting the 
hypothesis of self-organized criticality (SOC) in the solar activity, previously pro- 
posed for the solar flare occurrence, as discussed in the previous section. 

The observation that all the spectral peaks at periods of 2 years and shorter, 
here proposed to be harmonics of the QBO, have amplitudes decreasing nearly as 
1/f with frequency (clearly differently from the 11-year cycle, the harmonics of 
which decay much faster), also emerging in the same frequency range as the 1/f 
noise-like component (apart from the 2-year peak itself) indicates a connection of 
these oscillations to surface (photospheric) solar activity phenomena, such as the 
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Table 1: period (in years) of the spectral peaks in instantly maximal wavelet skeleton 
spectrum and those theoretically predicted from 11 and 2 year asymmetric oscilla- 
tors. 



Observed 


Theoretical 


28 ±2 




24 ±2 




19 ± 1 




16 ± 1 




13 ± 1 




11/n years 


10.8 ±0.8 


11 


8.8 ±0.7 




5.6 ±0.4 


5.5 


3.8 ±0.2 


3.67 


2.8 ±0.2 


2.75 


2/n years 


2.1 ±0.1 


2 


1.4 ±0.1 




1.07 ±0.08 


1 


0.97 ±0.06 




0.68 ±0.05 


0.67 


0.49 ± 0.04 


0.5 


0.42 ± 0.04 


0.4 


0.34 ±0.02 


0.33 




0.29 


0.26 ±0.02 


0.25 




0.22 


0.20 ±0.01 


0.2 
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solar flares. From theoretical point of view this observation again supports the 
model of double solar dynamo (Benevolenskaya, 1998, 2000) in which the 2-year 
oscillation is due to latitudal shear at the top of the convection zone, hence also the 
photosphere. The nature and properties of the surface dynamo producing both the 
QBO and intermittent bursts of l//-type (SOC), and the coupling with the main 
11-year dynamo at the basis of the convection zone is not within the scope of the 
present article, yet are to our opinion of great interest for solar physics and deserve 
further investigation. 

5 Conclusions 

We proposed the parallel use of two complementary skeleton spectra, extracted from 
the overall continuous wavelet spectrum, that are useful for studying the existence 
and evolution of stable or transient periodicities and noise-like components in time- 
series. As an indicative physical application, we analyzed the monthly averaged 
sunspot index. The results showed: 

(a) Clearer spectral peaks in the instantly maximal, compared to the overall, 
wavelet spectrum, that can be further attributed to the action of only two non- 
harmonic oscillators, namely the well-known 11-year and the 2-year (QBO), sup- 
porting the double-cycle solar dynamo model (Benevolenskaya, 1998, 2000). 

(b) The existence of a 1/ / noise background at timescales of 1 year and shorter, 
possibly being an imprint of the previously proposed self-organized critical state 
(SOC) of the solar magnetic fields and the intermittent energy releases, such as 
solar flares (e.g. Lu & Hamilton, 1991, Lu et al, 1993). 

We hope that the two wavelet skeleton spectra will be further useful in vari- 
ous signal processing applications, as data compression, filtering, data fitting and 
modeling. 
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Figure captions 

Figure 1. The monthly- averaged sunspot index, (b) The corresponding wavelet 
scalogram in grey color scale (whiter color corresponds to larger wavelet amplitudes), 
(c) The scale maximal and (d) the instantly maximal wavelet skeleton spectra. 

figure 2. The total wavelet power of the overall (upper curve), the instantly 
maximal (middle curve, the period in years of the spectral peaks is also noted) and 
scale maximal wavelet skeleton spectra in common, arbitrary units. 
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